function se = TUSIM_profse(b,Xc,h,X,use,resid)

n = size(X,1);
p = numel(b);
nuse = numel(use);
buse = b(use);
Xcu = Xc(:,use);
Xu = X(:,use);

prof_deriv = ((TUSIM_prof(buse*ones(1,nuse)+h*eye(nuse),Xcu) - TUSIM_prof(buse*ones(1,nuse)-h*eye(nuse),Xcu))/(2*h))';

S1 = (Xcu*buse).*(Xcu*buse > 0) - mean((Xcu*buse).*(Xcu*buse > 0));  % n x 1
Sb = Xu.*(resid*ones(1,size(Xu,2)))*inv(Xu'*Xu); %#ok<MINV>  n x pu
S2 = Sb*prof_deriv; % n x 1
S = S1/n+S2;

se = sqrt(S'*S);
